
%calibration
%estimating alp
cali=@(x)sum(grpsize.*( (ex_ante_allocation(x,updatek,ordering_list)-group).^2 ),'all'); 
x0=ones(size(group(1:end)))+0.5;
A = [];
b = [];
Aeq = [];
beq = [];
nonlcon=[];
lb= zeros(size(x0));
ub=[];
options = optimoptions(@fmincon,'Display','iter','Algorithm','sqp','MaxFunctionEvaluations',1e5,'MaxIterations',200,'OptimalityTolerance',1e-5,'StepTolerance',1e-4,'UseParallel',true);
%
tic
[alp,distance_pi]= fmincon(cali,x0,A,b,Aeq,beq,lb,ub,nonlcon,options)
toc
% 
alp_pi=alp;%[alp;5];
grpalp=alp;%[alp;5];
% % % % 
% % % calculating r^*
worstofffun=@(type)interim_total(alp_pi,type,ordering_list,grpsize)-1;
worstoff=fzero(worstofffun,[0.001,1]);
r_ast_pi=zeros(1,ngroup);
parfor agent=1:ngroup
     comb=ordering_list{agent}{1};
     lower=ordering_list{agent}{2};
     q_agent=ordering_list{agent}{3};
     numberpossibilities=ordering_list{agent}{4};
     r_ast_pi(agent)=interim_allocation(q_agent,grpalp,worstoff,comb,lower,numberpossibilities); 
end
r_ast_vec_pi=r_ast_pi(ikc); %1*570 
r_fac_pi=r_ast_vec_pi(ia2);
alp_fac_pi=alp_pi(ikc);
alp_fac_pi=alp_fac_pi(ia2);
r_ast_pi;
% % % 
% % % % obtaining the fitted expected allocation 
fitted_ea_pi=ex_ante_allocation(alp_pi,updatek,ordering_list);
ea_fac_vec_pi=fitted_ea_pi(ikc);
ea_fac_pi=ea_fac_vec_pi(ia2);
% % %
% % % revenue at r^*
r_ast_rev=revenue(alp_pi,updatek,ordering_list,worstoff*ones(ngroup,1),r_ast_pi,alp_fac_pi,grpsize)
% % % %revenue at r^sw
sw_wt=worstoff_types(updatek,alp_pi,ordering_list,fitted_ea_pi,alp_fac_pi)
r_sw_rev=revenue(alp_pi,updatek,ordering_list,sw_wt,fitted_ea_pi,alp_fac_pi,grpsize)
%propk_wt=worstoff_types(updatek,alp_pi,ordering_list,groupk/sum(fac_updatek),alp_fac_pi)
%r_propk_rev=revenue(alp_pi,updatek,ordering_list,propk_wt,groupk/sum(fac_updatek),alp_fac_pi,grpsize)
% % % % 
% % % % %%empirical revenue 
emp_wt=zeros(size(endowment));
emp_rev=zeros(nperiod,1);
parfor tt=1:nperiod
    empirical_share=min(max(endowment(tt,:),mini_demands'),fac_updatek'+mini_demands')-mini_demands';
    emp_wt(tt,:)=worstoff_types(updatek,alp_pi,ordering_list,empirical_share,alp_fac_pi)';
    emp_rev(tt)=revenue(alp_pi,updatek,ordering_list,emp_wt(tt,:)',empirical_share,alp_fac_pi,grpsize);
end
emp_rev

% test_share=zeros(1,57);
% test_share(find(fac_updatek==groupk(1),1,'first'))=1.7431-1.5691;
% %test_share(find(fac_updatek==groupk(4),2,'first'))=groupk(3);
% %test_share(find(fac_updatek==groupk(4),1,'last'))=1-grpsize(2)*groupk(2);
% test_wt=worstoff_types(updatek,alp_pi,ordering_list,test_share,alp_fac_pi)'
% test_rev=revenue(alp_pi,updatek,ordering_list,test_wt',test_share,alp_fac_pi,grpsize)
% 

% groupnumber=2;
% groupk(groupnumber)
% comb=ordering_list{groupnumber}{1};
% lower=ordering_list{groupnumber}{2};
% q_agent=ordering_list{groupnumber}{3};
% numberpossibilities=ordering_list{groupnumber}{4};
% interim_allocation(q_agent,alp_pi,0,comb,lower,numberpossibilities)
      
%%boundary solution: 
%inforents=info_rent(alp_pi,updatek,ordering_list)
corner_share=zeros(1,57);
corner_share(find(fac_updatek==groupk(1),1,'first'))=1;
%test_share(find(fac_updatek==groupk(4),2,'first'))=groupk(3);
%test_share(find(fac_updatek==groupk(4),1,'last'))=1-grpsize(2)*groupk(2);
corner_wt=worstoff_types(updatek,alp_pi,ordering_list,corner_share,alp_fac_pi)';
corner_rev=revenue(alp_pi,updatek,ordering_list,corner_wt',corner_share,alp_fac_pi,grpsize)

% % 
%trading volume, only account for the non-artifical facilities
cap_mean_eff_trading_vol=zeros(nperiod,1);
cap_eff_trading_vol_5=zeros(nperiod,1);
cap_eff_trading_vol_95=zeros(nperiod,1);
mean_eff_trading_vol=zeros(nperiod,1);
eff_trading_vol_5=zeros(nperiod,1);
eff_trading_vol_95=zeros(nperiod,1);
mean_ast=zeros(nperiod,1); ast_5=zeros(nperiod,1);ast_95=zeros(nperiod,1);
mean_sw=zeros(nperiod,1); sw_5=zeros(nperiod,1);sw_95=zeros(nperiod,1);
for t=1:nperiod
     alloc=endowment(t,:);%*supply_yearly.supply(t);
     adjalloc=min(max(endowment(t,:),mini_demands'),fac_updatek'+mini_demands');
     ast_alloc=r_fac_pi+mini_demands'; %1*57
     sw_alloc=ea_fac_pi'+mini_demands';
     ndraw=1000;
     nagent=nfac;
     vol=zeros(1,ndraw);
     voladj=zeros(1,ndraw);
     ast=zeros(1,ndraw);
     sw=zeros(1,ndraw);
     for rep=1:ndraw
        types=zeros(nagent,1);
        for agent=1:nagent
            a=alp_fac_pi(agent);
            types(agent)=invdist(rand(1),a); 
        end
        [~,ranking]=sort(types,'descend');
        allocation=realized_allocation(ranking,fac_updatek)'+mini_demands';
        vol(rep)=sum(abs(allocation-alloc))/2;%-abs(allocation(57)-alloc(57));
        voladj(rep)=sum(abs(allocation-adjalloc))/2;%-abs(allocation(57)-adjalloc(57));
        ast(rep)=sum(abs(allocation-ast_alloc))/2;%-abs(allocation(57)-ast_alloc(57));
        sw(rep)=sum(abs(allocation-sw_alloc))/2;%-abs(allocation(57)-sw_alloc(57));
     end
     mean_eff_trading_vol(t)=mean(vol);
     eff_trading_vol_5(t)=prctile(vol,5);
     eff_trading_vol_95(t)=prctile(vol,95);
     mean_ast(t)=mean(ast);
     ast_5(t)=prctile(ast,5); 
     ast_95(t)=prctile(ast,95);
     mean_sw(t)=mean(sw);
     sw_5(t)=prctile(sw,5); 
     sw_95(t)=prctile(sw,95);
     cap_mean_eff_trading_vol(t)=mean(voladj);
     cap_eff_trading_vol_5(t)=prctile(voladj,5);
     cap_eff_trading_vol_95(t)=prctile(voladj,95);
end


% % functions: 
% 

%v1
% function[meanvol,vol5,vol95,mean_ast,ast_5,ast_95]=trading_volume(k_agent,alloc,r,alp_agent) 
%     ndraw=1000;
%     nagent=length(k_agent);
%     vol=zeros(1,ndraw);
%     ast=zeros(1,ndraw);
%     alp=alp_agent;
%     for rep=1:ndraw
%         types=zeros(nagent,1);
%         for agent=1:nagent
%             a=alp(agent);
%             types(agent)=invdist(rand(1),a); 
%         end
%         [~,ranking]=sort(types,'descend');
%         allocation=realized_allocation(ranking,k_agent);
%         vol(rep)=sum(abs(allocation-alloc))/2;
%         ast(rep)=sum(abs(allocation-r))/2;
%     end
%     meanvol=mean(vol);
%     vol5=prctile(vol,5);
%     vol95=prctile(vol,95);
%     mean_ast=mean(ast);ast_5=prctile(ast,5);ast_95=prctile(ast,95);
% end
function[rev]=revenue(grpalp,grpk,ordering_list,worstofftypes,alloc,alpfac,grpsize)
     nagent=length(alloc);
     ngroup=length(grpk);
     rev=zeros(nagent,1);
     for agent=1:nagent
         if ngroup==nagent
            groupnumber=agent;
        else
            groupnumber=find(grpalp==alpfac(agent));
        end
        comb=ordering_list{groupnumber}{1};
        lower=ordering_list{groupnumber}{2};
        q_agent=ordering_list{groupnumber}{3};
        numberpossibilities=ordering_list{groupnumber}{4};
        a=grpalp(groupnumber);
        w=worstofftypes(agent);
        fun2=@(type)interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities).*buyertypetimesdensity(a,type);
        fun1=@(type)interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities).*sellertypetimesdensity(a,type);
        rev(agent)=integral(fun1,0,w)+integral(fun2,w,1)-max(0,min(alloc(agent),grpk(groupnumber)))*w;
%         m=10;
%         x0=0;
%         x1=1;
%         dx=(x1-x0)/m;
%         int=0;
%         parfor i=1:m
%             int=int+ fun(x0+(i-0.5)*dx)*dx
%         end
%         eq(agent)=int;
     end
     if ngroup==nagent
        rev=sum(rev.*grpsize);
     else
        rev=sum(rev);
     end
end

function[wt]=worstoff_types(groupk,grpalp,ordering_list,alloc,alpfac) % note alloc has to be the same size as grpk
    nagent=length(alloc);
    ngroup=length(groupk);
    wt=zeros(nagent,1);
    for agent=1:nagent
        if ngroup==nagent
            groupnumber=agent;
        else
            groupnumber=find(grpalp==alpfac(agent));
        end
        comb=ordering_list{groupnumber}{1};
        lower=ordering_list{groupnumber}{2};
        q_agent=ordering_list{groupnumber}{3};
        numberpossibilities=ordering_list{groupnumber}{4};
        fun=@(type)min(alloc(agent),groupk(groupnumber))-interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
        if fun(0)<0
             wt(agent)=0;
        elseif fun(1)>0
            wt(agent)=1;
        else
            wt(agent)=fzero(fun,[0,1]);
        end
    end
end

function[eq]=ex_ante_allocation(grpalp,grpk,ordering_list)
    nagent=length(grpk);
    eq=zeros(nagent,1);
    for agent=1:nagent
        a=grpalp(agent);
        comb=ordering_list{agent}{1};
        lower=ordering_list{agent}{2};
        q_agent=ordering_list{agent}{3};
        numberpossibilities=ordering_list{agent}{4};
        fun=@(type)density(a,type).*interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
        eq(agent)=integral(fun,0,1);
%         m=10;
%         x0=0;
%         x1=1;
%         dx=(x1-x0)/m;
%         int=0;
%         parfor i=1:m
%             int=int+ fun(x0+(i-0.5)*dx)*dx
%         end
%         eq(agent)=int;
    end
end
function[eq]=info_rent(grpalp,grpk,ordering_list)
    nagent=length(grpk);
    eq=zeros(nagent,1);
    for agent=1:nagent
        a=grpalp(agent);
        comb=ordering_list{agent}{1};
        lower=ordering_list{agent}{2};
        q_agent=ordering_list{agent}{3};
        numberpossibilities=ordering_list{agent}{4};
        fun=@(type)interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
        eq(agent)=integral(fun,0,1);
%         m=10;
%         x0=0;
%         x1=1;
%         dx=(x1-x0)/m;
%         int=0;
%         parfor i=1:m
%             int=int+ fun(x0+(i-0.5)*dx)*dx
%         end
%         eq(agent)=int;
    end
end

function[suma]=interim_total(grpalp,type,ordering_list,grpsize) %remember to change interim total 
    nagent=length(grpalp);
    suma=0;
    parfor j=1:nagent
        comb=ordering_list{j}{1};
        lower=ordering_list{j}{2};
        q_agent=ordering_list{j}{3};
        numberpossibilities=ordering_list{j}{4};
        suma=suma+grpsize(j)*interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities);
    end    
end

%v1:
% function[agent]=interim_allocation(grp,grpalp,type,grpk,comb,lower) %alp needs to be a col vector   
%     q_agent=agent_allocation(grp,comb,grpk);
%     alp=grpalp;
%     higher=comb;
%     ntype=length(type);
%     agent=zeros(size(type));
%     parfor tt=1:ntype
%         prob11=prob1(alp,type(tt),higher);
%         prob22=prob2(alp,type(tt),lower);
%         agent(tt)=sum(prob11.*prob22.*q_agent,'all');
%     end
% end

%v2: use v2 of ordering_list
function[agent]=interim_allocation(q_agent,grpalp,type,comb,lower,numberpossibilities) %alp needs to be a col vector   
    %q_agent=agent_allocation(grp,comb,grpk);
    alp=grpalp;
    higher=comb;
    ntype=length(type);
    agent=zeros(size(type));
    parfor tt=1:ntype
        prob11=prob1(alp,type(tt),higher);
        prob22=prob2(alp,type(tt),lower);
        agent(tt)=sum(prob11.*prob22.*q_agent.*numberpossibilities,'all');
    end
end

function[q]=realized_allocation(order,fac_updatek)
    nagent=length(order);
    q=zeros(nagent,1);
    rest=1;
    allocated=0;
    k=fac_updatek;
    for rk=1:nagent
        agent=order(rk);
        q(agent)=min(rest,k(agent));
        allocated=allocated+q(agent);
        rest=max(0,1-allocated);
    end
end


function[z1]=density(a,type)
    z1= a.*(type).^(a-1);
end

function[z1]=buyertypetimesdensity(a,type)
    z1= a.*(type).^(a)-(1-type.^a);
end

function[z1]=sellertypetimesdensity(a,type)
    z1= a.*(type).^(a)+(type.^a);
end

function[t]=invdist(u,a) %given quantile, output type
    t= u.^(1./a);
end
function[z2]=prob1(a,type,higher)
    z2=prod((1-(type).^a).^higher,1)';
end

function[z3]=prob2(a,type,lower)
    z3=prod( ( (type).^a ).^lower,1 )';
end


